home *** CD-ROM | disk | FTP | other *** search
/ Our Solar System / Our Solar System.iso / miscprog / librat / oprecess.c < prev    next >
C/C++ Source or Header  |  1990-07-31  |  9KB  |  419 lines
  1. /* opreces.c
  2.  * Rigorous precession of orbital elements
  3.  *
  4.  * See euler.txt for explanations.
  5.  * Program by Steve Moshier.
  6.  */
  7. #define SALONE 0
  8. double sin(), cos(), zatan2(), fabs();
  9. #define N 3
  10.  
  11.  
  12. #if SALONE
  13.  
  14. double J2000 = 2451545.0; /* 2000 January 1.5 */
  15. double STR = 4.8481368110953599359e-6; /* radians per arc second */
  16. double DTR = 1.7453292519943295769e-2; /* radians per degree */
  17. double RTD = 5.7295779513082320877e1; /* degrees per radian */
  18. double PI = 3.14159265358979323846;
  19.  
  20. /* Indicate whether input is argument of perihelion (1)
  21.  * or longitude of perihelion (0):
  22.  */
  23. #define ARGPERIH 1
  24.  
  25. main()
  26. {
  27. double JD1, JD2, e[3];
  28.  
  29. printf( "Enter Julian date of given orbital elements ? " );
  30. scanf( "%lf", &JD1 );
  31. printf( "%.2lf\n", JD1 );
  32. printf( "Julian date of desired orbital elements ? " );
  33. scanf( "%lf", &JD2 );
  34. printf( "%.2lf\n", JD2 );
  35.  
  36. /* Find the rotation from the ecliptic of JD1 to the orbit.
  37.  */
  38. printf( "Enter longitude of ascending node, in degrees ? " );
  39. scanf( "%lf", &e[0] );
  40. printf( "%.6lf\n", e[0] );
  41. printf( "Inclination ? " );
  42. scanf( "%lf", &e[1] );
  43. printf( "%.6lf\n", e[1] );
  44.  
  45. #if ARGPERIH
  46. printf( "Arg perihelion ? " );
  47. scanf( "%lf", &e[2] );
  48. #else
  49. printf( "Longitude of the perihelion ? " );
  50. scanf( "%lf", &e[2] );
  51. printf( "%.6lf\n", e[2] );
  52. e[2] -= e[0];
  53. #endif
  54.  
  55. e[2] *= DTR; /* convert degrees to radians */
  56. e[1] *= DTR;
  57. e[0] *= DTR;
  58.  
  59. oprecess( JD1, JD2, e );
  60.  
  61. printf( "\nElements for equinox of %.2lf\n", JD2 );
  62. printf( "Node = %.6lf\n", RTD*e[0] );
  63. printf( "Inclination = %.6lf\n", RTD*e[1] );
  64. printf( "Arg perihelion = %.6lf\n", RTD*e[2] );
  65. printf( "Longitude of the perihelion = %.6lf\n", RTD*(e[2]+e[0]) );
  66. }
  67.  
  68.  
  69.  
  70. /* Matrix transpose.
  71.  * B, the output, can occupy the same storage locations as
  72.  * A, the input.
  73.  */
  74. mtransp( A, B )
  75. double A[], B[];
  76. {
  77. double x, y;
  78. int r, c, rc, cr;
  79.  
  80. for( r=0; r<N; r++ )
  81.     {
  82.     for( c=0; c<=r; c++ )
  83.         {
  84.         rc = N*r+c;
  85.         cr = N*c+r;
  86.         x = A[rc];
  87.         y = A[cr];
  88.         B[rc] = y;
  89.         B[cr] = x;
  90.         }
  91.     }
  92. }
  93.  
  94.  
  95.  
  96.  
  97. /* 3x3 matrix multiply
  98.  * C = A B
  99.  * C may occupy the same storage as either A or B.
  100.  */
  101. mmpy3( A, B, C )
  102. double A[3][3], B[3][3], C[3][3];
  103. {
  104. double ans[3][3];
  105. double s;
  106. int i, r, c, k;
  107.  
  108. for( r=0; r<3; r++ )
  109.     {
  110.     for( c=0; c<3; c++ )
  111.         {
  112.         s = 0.0;
  113.         for ( i=0; i<3; i++ )
  114.             s += A[r][i] * B[i][c];
  115.         ans[r][c] = s;
  116.         }
  117.     }
  118.  
  119. for( r=0; r<3; r++ )
  120.     {
  121.     for( c=0; c<3; c++ )
  122.         {
  123.         C[r][c] = ans[r][c];
  124.         }
  125.     }
  126. }
  127.  
  128.  
  129.  
  130. /* Construct the matrix that represents the composite of
  131.  * successive rotations by the three Euler angles:
  132.  *   first by an angle phi about the z axis,
  133.  *   second by an angle theta about the new x axis,
  134.  *   third by an angle psi about the final z axis.
  135.  * The input array e[] contains the three angles in the order
  136.  * phi, theta, psi.  The 3x3 output matrix M[][] is the 
  137.  * composite transformation of the three rotations in
  138.  * Cartesian coordinates.
  139.  *   The matrix elements have been written out directly
  140.  * in trigonometric form so that the derivation of the
  141.  * inverse function MtoE() will be clear.
  142.  */
  143. EtoM( e, M )
  144. double e[3];
  145. double M[3][3];
  146. {
  147. double a, b;
  148. double sinpsi, cospsi, sinth, costh, sinphi, cosphi;
  149.  
  150. sinpsi = sin(e[2]);
  151. cospsi = cos(e[2]);
  152. sinth = sin(e[1]);
  153. costh = cos(e[1]);
  154. sinphi = sin(e[0]);
  155. cosphi = cos(e[0]);
  156. a = costh*sinphi;
  157. b = costh*cosphi;
  158. M[0][0] = cospsi*cosphi - a*sinpsi;
  159. M[0][1] = cospsi*sinphi + b*sinpsi;
  160. M[0][2] = sinpsi*sinth;
  161. M[1][0] = -sinpsi*cosphi - a*cospsi;
  162. M[1][1] = -sinpsi*sinphi + b*cospsi;
  163. M[1][2] = cospsi*sinth;
  164. M[2][0] = sinth*sinphi;
  165. M[2][1] = -sinth*cosphi;
  166. M[2][2] = costh;
  167. }
  168.  
  169.  
  170.  
  171.  
  172. /* Deduce the three Euler angles e[]
  173.  * from the input coordinate rotation matrix M[][].
  174.  * This is done by trigonometry from the elements of M.
  175.  *
  176.  * Since
  177.  *   M[2][1] = -sin(theta) * cos(phi)
  178.  *   M[2][0] =  sin(theta) * sin(phi),
  179.  * then phi = arctan( M[2][0]/M[2][1].
  180.  *
  181.  * Similarly,
  182.  *  M[0][2] = sin(psi) * sin(theta)
  183.  *  M[1][2] = cos(psi) * sin(theta),
  184.  * so psi = arctan( M[0][2]/M[1][2] ).
  185.  *
  186.  * Finally,
  187.  *   M[2][2] = cos(theta)
  188.  *   M[1][2]/cos(psi) = sin(theta),
  189.  * which gives the arctangent of theta.
  190.  */
  191. MtoE( M, e )
  192. double M[3][3], e[];
  193. {
  194. double a, b, phi, theta, psi;
  195. double zatan2(), fabs(), sqrt(), sin(), cos();
  196.  
  197. phi = zatan2( -M[2][1], M[2][0] );
  198. a = M[0][2];
  199. b = M[1][2];
  200. psi = zatan2( b, a );
  201.  
  202. if( fabs(a) > fabs(b) )
  203.     a = a/sin(psi);
  204. else
  205.     a = b/cos(psi);
  206.  
  207. theta = zatan2( M[2][2], a );
  208.  
  209. e[0] = phi;
  210. e[1] = theta;
  211. e[2] = psi;
  212. }
  213.  
  214.  
  215.  
  216.  
  217.  
  218.  
  219. /*                            zatan2()
  220.  *
  221.  *    Quadrant correct inverse circular tangent
  222.  *
  223.  *
  224.  *
  225.  * SYNOPSIS:
  226.  *
  227.  * double x, y, z, zatan2();
  228.  *
  229.  * z = zatan2( x, y );
  230.  *
  231.  *
  232.  *
  233.  * DESCRIPTION:
  234.  *
  235.  * Returns radian angle between 0 and +2pi whose tangent
  236.  * is y/x.
  237.  *
  238.  *
  239.  *
  240.  * ACCURACY:
  241.  *
  242.  * See atan.c.
  243.  *
  244.  */
  245.  
  246. /*
  247. Cephes Math Library Release 2.0:  April, 1987
  248. Copyright 1984, 1987 by Stephen L. Moshier
  249. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  250. Certain routines from the Library, including this one, may
  251. be used and distributed freely provided this notice is retained
  252. and source code is included with all distributions.
  253. */
  254.  
  255. /* include "prec.h" */
  256.  
  257.  
  258. double zatan2( x, y )
  259. double x, y;
  260. {
  261. double z, w;
  262. short code;
  263. double atan();
  264.  
  265.  
  266. code = 0;
  267.  
  268. if( x < 0.0 )
  269.     code = 2;
  270. if( y < 0.0 )
  271.     code |= 1;
  272.  
  273. if( x == 0.0 )
  274.     {
  275.     if( code & 1 )
  276.         return( 1.5*PI );
  277.     if( y == 0.0 )
  278.         return( 0.0 );
  279.     return( 0.5*PI );
  280.     }
  281.  
  282. if( y == 0.0 )
  283.     {
  284.     if( code & 2 )
  285.         return( PI );
  286.     return( 0.0 );
  287.     }
  288.  
  289.  
  290. switch( code )
  291.     {
  292.     case 0: w = 0.0; break;
  293.     case 1: w = 2.0 * PI; break;
  294.     case 2:
  295.     case 3: w = PI; break;
  296.     }
  297.  
  298. z = atan( y/x );
  299.  
  300. return( w + z );
  301. }
  302.  
  303. #else
  304. extern double J2000, STR, DTR, RTD, PI;
  305. #endif /* SALONE */
  306.  
  307.  
  308. /* Subroutine to precess orbital elements
  309.  * from Julian date JD1 to Julian date JD2.
  310.  * The input elements, corresponding to precessional epoch JD1, are
  311.  * e[0] = node
  312.  * e[1] = inclination
  313.  * e[2] = arg perihelion
  314.  */
  315. oprecess( JD1, JD2, e )
  316. double JD1, JD2, e[];
  317. {
  318. double x, max;
  319. double P[3][3], Q[3][3];
  320.  
  321. /* Find the precession matrix to go from JD2 to JD1.
  322.  */
  323. epmat( JD1, P ); /* J2000 to JD1 */
  324. epmat( JD2, Q ); /* J2000 to JD2 */
  325. mtransp( Q, Q ); /* transpose = JD2 to J2000 */
  326. mmpy3( P, Q, P ); /* JD2 to J2000 to JD1 */
  327.  
  328. EtoM( e, Q ); /* Convert input Euler angles to rotation matrix */
  329.  
  330. /* Calculate the composite rotation
  331.  * ecliptic JD2 -> ecliptic JD1 -> orbit.
  332.  */
  333. mmpy3( Q, P, Q );
  334.  
  335. /* Find the Euler angles of the result
  336.  */
  337. MtoE( Q, e );
  338. }
  339.  
  340. /* The following subprogram, epmat( JD, P ),
  341.  * constructs the precession matrix in ecliptic coordinates
  342.  * to go from the epoch J2000.0 to the date JD.  Output P is the
  343.  * 3x3 rotation matrix.  To go in the opposite direction,
  344.  * from JD to J2000.0, the required matrix is just the transpose of P.
  345.  * Note that the obliquity of date is not required here since the
  346.  * coordinates are ecliptic, not equatorial.
  347.  */
  348.  
  349. /* Accumulated precession in longitude.  Coefficients are from:
  350.  * J. Laskar, "Secular terms of classical planetary theories
  351.  * using the results of general theory," Astronomy and Astrophysics
  352.  * 157, 59070 (1986).
  353.  */
  354. static double pAcof[] = {
  355.  -8.66e-10, -4.759e-8, 2.424e-7, 1.3095e-5, 1.7451e-4, -1.8055e-3,
  356.  -0.235316, 0.07732, 111.1971, 50290.966 };
  357.  
  358. /* Node and inclination of the earth's orbit computed from
  359.  * Laskar's data as explained in the following paper:
  360.  * P. Bretagnon and G. Francou, "Planetary theories in rectangular
  361.  * and spherical variables. VSOP87 solutions," Astronomy and
  362.  * Astrophysics 202, 309-315 (1988).
  363.  */
  364. static double nodecof[] = {
  365. 6.6402e-16, -2.69151e-15, -1.547021e-12, 7.521313e-12, 6.3190131e-10, 
  366. -3.48388152e-9, -1.813065896e-7, 2.75036225e-8, 7.4394531426e-5,
  367. -0.042078604317, 3.052112654975 };
  368.  
  369. static double inclcof[] = {
  370. 1.2147e-16, 7.3759e-17, -8.26287e-14, 2.503410e-13, 2.4650839e-11, 
  371. -5.4000441e-11, 1.32115526e-9, -5.998737027e-7, -1.6242797091e-5,
  372.  0.002278495537, 0.0 };
  373.  
  374.  
  375. epmat( JD, P )
  376. double JD;
  377. double P[];
  378. {
  379. double T, s, pA, i, Omega;
  380. double e[3];
  381. double *p;
  382. int j;
  383.  
  384. T = (JD - J2000)/365250.0; /* thousands of years from J2000.0 */
  385.  
  386. /* Accumulated precession in longitude
  387.  */
  388. p = pAcof;
  389. pA = *p++;
  390. for( j=0; j<9; j++ )
  391.     pA = pA * T + *p++;
  392. pA *= STR * T;
  393.  
  394. /* Node of the moving ecliptic on the J2000 ecliptic.
  395.  */
  396. p = nodecof;
  397. Omega = *p++;
  398. for( j=0; j<10; j++ )
  399.     Omega = Omega * T + *p++;
  400.  
  401. /* Inclination of the moving ecliptic to the J2000 ecliptic.
  402.  */
  403. p = inclcof;
  404. i = *p++;
  405. for( j=0; j<10; j++ )
  406.     i = i * T + *p++;
  407.  
  408. /* Set up the Euler angles */
  409. e[0] = Omega;
  410. e[1] = i;
  411. e[2] = -(Omega + pA );
  412.  
  413. /* Construct the rotation matrix
  414.  * to go from J2000 to JD
  415.  */
  416. EtoM( e, P );
  417. }
  418.  
  419.